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Abstract. We describe the use of pyramid decomposition in Normaliz, a software tool 
for the computation of Hilbert bases and enumerative data of rational cones and affine 
monoids. Pyramid decomposition in connection with efficient parallelization and stream- 
lined evaluation of simplicial cones has enabled Normaliz to process triangulations of size 
« 5- 10 11 that arise in the computation of Hilbert series related to combinatorial voting 
theory. 



1. Introduction 

Normaliz is a software tool for the computation of Hilbert bases and enumerative data 
of rational cones and affine monoids. In the 14 years of its existence it has found numerous 
applications; for example, see Bogart, Raymond and Thomas 0, Craw, Maclagan and 
Thomas 0H, Kappl, Ratz and Staudt EB or Sturmfels and Welker (22). Normaliz is 
used in polymake lfT6ll and Regina O. 

The mathematics of the very first version was described in [8 1, and the ideas leading to 
version 2.2 (2009) are contained in [|6). In this article we want to document the most recent 
development that has extended the scope of Normaliz by several orders of magnitude. 

From the very beginning Normaliz has used lexicographic triangulations. (It also con- 
tains a triangulation free Hilbert basis algorithm; see [|6]). Lexicographic triangulations 
are essentially characterized by being incremental in the following sense. Suppose that the 

cone C is generated by vectors xi , . . . ,x n 6 Mr and set C, = JLpci H hM + jc,-, i = 0, . . . ,n. 

Then the lexicographic triangulation A (for the ordered system x\ , . . . ,x n ) restricts to a tri- 
angulation of C, for i = 0, ... ,n. Lexicographic triangulations are easy to compute, and 
go very well with Fourier-Motzkin elimination that computes the support hyperplanes of 
C by successive extension from Q to Q+i, i = 0,...,n— 1. The triangulation A, of Q is 
extended to C i+ i by all simplicial cones F + HLpCj+i where F 6 A; is visible from Xi + \. 

As simple as the computation of the lexicographic triangulation is, the algorithm in 
the naive form just described has two related drawbacks: (i) one must store A, and this 
becomes very difficult for sizes > 10 8 ; (ii) in order to find the facets F that are visible from 
Xj + { we must match the simplicial cones in A ; with the support hyperplanes of C, that are 
visible from jc,-+i. While (i) is a pure memory problem, (ii) quickly leads to impossible 
computation times. 

Pyramid decomposition is the basic idea that has enabled Normaliz to compute dimen- 
sion 24 triangulations of size ss 5 ■ 10 11 in acceptable time on standard multiprocessor 
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systems such as SUN xFire 4450 or Dell PowerEdge R910. Instead of going for the lex- 
icographic triangulation directly, we first decompose C into the pyramids generated by 
xi + \ and the facets of Q that are visible from i = 0, . . . ,n— 1. These pyramids (of 
level 0) are then decomposed into pyramids of level 1 etc. While the level decomposi- 
tion need not be a polyhedral subdivision in the strict sense, pyramid decomposition stops 
after finitely many iterations at the lexicographic triangulation (see Section 12]). 

Pure pyramid decomposition is extremely memory friendly, but its computation times 
are even more forbidding than those of pure lexicographic triangulation since too many 
Fourier-Motzkin eliminations become necessary, and almost all of them are inevitably 
wasted. That Normaliz can nevertheless cope with extremely large triangulations relies 
on a well balanced combination of both strategies that we outline in Section [3] 

In Section |4] we describe the steps by which Normaliz evaluates the simplicial cones in 
the triangulation for the computation of Hilbert bases, volumes and Hilbert series. The 
evaluation almost always takes significantly more time than the triangulation. Therefore 
it must be streamlined as much as possible. For the Hilbert series Normaliz uses a Stanley 
decomposition ETI . That it can be found efficiently relies crucially on an idea of Koppe 
and Verdoolaege [18J. 

We document the scope of Normaliz' computations in Section |5J Our main examples 
come from combinatorial voting theory that we found in Schurmann's paper [20] . The de- 
sire to master the Hilbert series computations asked for in fl20l was an important stimulus 
in the recent development of Normaliz. 

For Hilbert basis computations pyramid decomposition has a further and sometimes 
tremendous advantage: one can avoid the triangulation of those pyramids for which it is 
a priori clear that they will not supply new candidates for the Hilbert basis. This observa- 
tion, on which the contribution of the authors to [5] is based, triggered the use of pyramid 
decomposition as a general principle. See Remark[8]for a brief discussion and Section [531 
for data of computations. 

It is an important aspect of pyramid decomposition that it is very parallelization friendly 
since the pyramids can be treated independently of each other. Normaliz uses OpenMP 
for shared memory systems. Needless to say that triangulations of the size mentioned 
above can hardly be reached in serial computation. 

2. Lexicographic triangulation and Pyramid decomposition 

Consider vectors xi,...,x n G M. d . For Normaliz these must be integral vectors, but 
integrality is irrelevant in this section. We want to compute a triangulation of the cone 

C = cone(jq , . . . ,x n ) = M. + xi H h M. + x n 

with rays through x\,... ,x n . Such a triangulation is a polyhedral subdivision of C into 
simplicial subcones o generated by linearly independent subsets of {x\,...,x n }. 
For a triangulation E of a cone C and a subcone C' we set 

I\Cf = {anCf:aeZ, dim a DC' = dimC'}. 

In general E|C ; need not be a triangulation of C', but it is so if C' is a face of C. 

The lexicographic (or placing) triangulation A(xi,. . . ,x n ) of cone(jci, . . . ,x n ) can be 
defined recursively as follows: (i) the triangulation of the zero cone is the trivial one, 
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(ii) A(x\ , . . . ,x n ) is given by 

A(xi,...,x„) = A(xi, . . . ,x M _i) U {cone(a,^„) : o e A(xi, . . . ,x n _i) visible from x n } 

where c is visible from x„ if x n £ cone(^i, . . . ,x n _i) and the line segment [x n ,y] for every 
point y of a intersects cone(^i, . . . ,x n _i) only in y. Note that a polyhedral complex is 
always closed under the passage to faces, and the definition above takes care of it. In 
the algorithms below a polyhedral subdivision can always be represented by its maximal 
faces which for convex full dimensional polyhedra are the full dimensional cones in the 
subdivision. For simplicial subdivisions of cones one uses of course that the face structure 
is completely determined by set theory: every subset E of the set of generators spans a 
conical face of dimension 

We state some useful properties of lexicographic triangulations: 

Proposition 1. With the notation introduced, let Ci = cone(xi ,...,.*,•) and A,- = A(xi , . . . , 
Xi )fori = 

(1) A n is the unique triangulation ofC with rays through xi, . . . ,x n that restricts to a 
triangulation ofCtfor i = 1 , . . . , n and satisfies A\Q = A|C,-_i ifQ = C,-_i. 

(2) For every face F ofC the restriction A\F is the lexicographic triangulation A(xt 1 , 
. . . ,x im ) where {x h , . . . ,x im } = FD {xi, . . . ,x n } and h < • • • < i m - 

(3) T/'dimQ > dimQ i, then A = A(xi, . . . ,x;_2,x;,x;-i,xh-i, . . . ,x n ). 

(4) A = A(xjj, . . . ,Xi d ,Xj r , . . . ,Xj nd ) where (ii,...,id) is the lexicographic smallest 
index vector of a rank d subset of {xi, . . . ,x n } and j\ < ■■ ■ < j n -d ^ sts tne com- 
plementary indices. 

Proof (1) By construction it is clear that A n satisfies the properties of which we claim that 
they determine A uniquely. On the other hand, the extension of A,_i to a triangulation of 
Q is uniquely determined if one does not introduce further rays: the triangulation of the 
part V of the boundary of Q _i that is visible from xi has to coincide with the restriction 
of Aj-i toV. 

(2) One easily checks that A|F satisfies the conditions in (1) that characterize A(jc ; , , . . . , 

X im)- 

(3) It is enough to check the claim for i = n. Then the only critical point for the 
conditions in (1) is whether A(xi, . . . ,x„_2,x n ,x n _i) restricts to C n -\. But this is the 
case since C„_i is a facet of C if dimC > dimC„_i. 

(4) follows by repeated application of (3). □ 

In the following we will assume that C is full dimensional: dimC = d = &\raR d . Part 
(4) helps us to keep the data structure of lexicographic triangulations simple: right from 
the start we need only to work with the list of dimension d simplicial cones of A by search- 
ing x/j , . . . ,Xi d first, choosing cone(jc ;i , . . . ,Xi d ) as the first J-dimensional simplicial cone 
and subsequently extending the list as prescribed by the definition of the lexicographic 
triangulation. In other words, we can assume that x\, . . . are linearly independent, and 
henceforth we will do so. 

In order to extend the triangulation we must of course know which facets of Q_i are 
visible from x,-. Recall that a cone C of dimension d in R d has a unique irredundant 
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representation as an intersection of linear half spaces: 

where Jf(C) is a finite set of hyperplanes and the orientation of the closed half spaces 
H and H+ is chosen in such a way that C C H + for H G Jff(C). For H G Jf(Q-i) the 
facet H fl Ci-\ is visible from Xj if and only if x\ lies in the open half space H < = H \ H. 
When we refer to support hyperplane in the following we always mean those that appear 
in the irredundant decomposition of C since only they are important in the algorithmic 
context. 

Hyperplanes are represented by linear forms X G (IR rf )*, and we always work with the 
basis e\,...,e^ that is dual to the basis e\, . . . , of unit vectors. For rational hyperplanes 
the linear form X can always be chosen in such a way that it has integral coprime co- 
efficients and satisfies X (x) > for x G C. This choice determines X uniquely. (If one 
identifies e*,...,e*j with ei, . . . , ej via the standard scalar product, then X is nothing but 
the primitive integral inner (with respect to C) normal vector of H.) For later use we 
define the (lattice) height of x G E over// by 

ht H (x) = |A(x)|. 

If F = CDH is the facet of C cut out by H, we set htp(x) = ht H (x). 

We can now describe the computation of the triangulation A(xi, . . . ,x n ) in a more for- 
mal way in Tabled] For simplicity we will identify a simplicial cone o with its generating 
set C {xi, .. . ,x n }. It should be clear from the context what is meant. For further use we 
introduce the notation 

J^*(C,x) = {H G Jf(C),x G H*} where * G {<,>,+,-}. 

Table[T] formalizes the computation of the lexicographic triangulation. 

LexTriangulation (xi , . . . , x n ) ExtendTri (/) 

1 ADDS implex (xi, . . . ,Xd) 1 parallel for H G Jt < (Ci-i,Xi) 

2 for i 4- d + 1 to n 2 do 

3 do 3 for <7 G A(xi,...,x;_i) 

4 ExtendTri (?) 4 do 

5 FindNewHyp(i) 5 if |<rDH| =</- 1 

6 then AddSimplex(xjU (aHH)) 

Table 1 . Incremental building of lexicographic triangulation 



The function AddS implex adds a simplicial cone to the (initially empty) list of sim- 
plicial cones that, upon completion, contains the lexicographic triangulation of C. The 
function FindNewHyp computes J$?(Ci) from Jt?(Ci-i) by Fourier-Motzkin elimina- 
tion. (It does nothing if xi G C,-_i.) Its Normaliz implementation has been described 
in great detail in [6]; therefore we skip it here. The function ExtendTri does exactly 
what its name says : it extends the triangulation A(xi , . . . , x ( -_ i ) of Q_ i to the triangulation 
A(xi, . . . ,Xj) of Q (again doing nothing if x/ G Q_i). 
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Note that the set of hyperplanes over which the loop in ExtendTri runs is given by 

Jf<{Q- h Xi). 

One is tempted to improve ExtendTri by better bookkeeping and using extra infor- 
mation on triangulations of cones. We discuss our more or less fruitless attempts in the 
following remark. 

Remark 2. (a) If one knows the restriction of A(xi, . . . ,X(_i) to the facets of C,-_i, then 
A(xi, . . . ,X() can be computed very fast. However, unless i = n, the facet triangulation 
must now be extended to the facets of Q, and this step eats up the previous gain, as exper- 
iments have shown, at least for the relatively small triangulations to which ExtendTri 
is really applied after the pyramid decomposition described below. 

(b) The test of the condition \o DH\ = d — 1 is satisfied if and only if d — 1 of the 
generators of o lie in H. Its verification can be accelerated if one knows which facets of 
the J-dimensional cones in A(xi, . . . ,x^_i) are already shared by another simplicial cone 
in A(xi, . . . ,x,'_i), and are therefore not available for the formation of a new simplicial 
cone. But the extra bookkeeping requires more time than is gained by its use. 

(c) One refinement is used in our implementation, though its influence is almost unmea- 
surable. Each simplicial cone in A(xi, . . . ,x^_i) has been added with a certain generator 
Xj, j < i. (The first cone is considered to be added with each of its generators.) It is 
not hard to see that only those simplicial cones that have been added with a generator 
xj E H can satisfy the condition |oTW| = d — 1, and this information is used to reduce 
the number of pairs (//, a) to be tested. 

(d) If \HCl {xi, . .. ,Xj_i}| = d — 1, then H 6 J^ < (C i _i,x i ) produces exactly one new 
simplicial cone of dimension d, namely cone(x;,// fl {xi, . . . ,x,-_i}), and therefore the 
loop over o can be suppressed. 

The product |J^ < (Q_i,x i -)| • |E| determines the complexity of ExtendTri. Even 
though the loop over H is parallelized (as indicated by parallel for), the time spent in 
ExtendTri can be extremely large. (The "exterior" loops in FindNewHyp are par- 
allelized as well.) The second limiting factor for ExtendTri is memory: it is already 
critical to store triangulations of size 10 8 and impossible for size > 10 9 . Therefore the 
direct approach to lexicographic triangulations does not work for truly large cones. 

Now we present a radically different way to lexicographic triangulations via iterated 
pyramid decompositions. The cones that appear in this type of decomposition are called 
pyramids since their cross-section polytopes are pyramids in the usual sense, namely of 
type conv(F,x) where F is a facet and x is a vertex not contained in F. 

Definition 3. The pyramid decomposition n(xi, . . . ,x M ) of C = cone(;q, . . . ,x n ) is recur- 
sively defined as follows: it is the trivial decomposition for n = 0, and 

n(xi,...,x n ) =n(jci,...,jc„_i)U{cone(i 7 ,Jc„) : 

F a face of C(x\, . . . ,x n -i) visible from x„}. 

Note that the pyramid decomposition is not a polyhedral subdivision in the strong sense: 
the intersection of two faces F and F' need not be a common face of F and F' (but is 
always a face of F or F'). See Figure[T]for an example. 
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Figure 1 . Cross-section of a pyramid decomposition 



In order to iterate the pyramid decomposition we set il°(jci, . . . ,x n ) = Tl(x\, . . .,x n ), 
and 

n\x h ...,x n )= |J {nixr.x.eP)} forfc>o. 

Pen*- 1 

Note that this recursion cannot descend indefinitely, since the number of generators goes 
down in each recursion level. We denote the total pyramid decomposition by n°°(xi , . . . , 
x n ). More precisely: 

Proposition 4. One has TY°(x\, . . . ,x n ) = Tl n ~ d (x\, . ..,x n )= A(jq, . . . ,x n ). 

Proof. In the case n = d, the pyramid decomposition is obviously the face lattice of C, 
and therefore coincides with the lexicographic triangulation. For n > d the first full di- 
mensional pyramid reached is the simplicial cone cone (x\ , . . . , Xd) ■ All the other pyramids 
have at most n — 1 generators, and so we can use induction: For each P 6 II (*i, . . . , x n ) 
the total pyramid decomposition of P is the lexicographic triangulation A(x, : x\ G P). 
According to Proposition [TJ2) these triangulations match along the common boundaries 
of the pyramids, and therefore constitute a triangulation of C. It evidently satisfies the 
conditions in Proposition [T]T). □ 

This leads to a recursive computation of A(^i , . . . ,x n ) by the algorithms in Table |2l The 

TotalPyrDec(xi, . . . ,x n ) ProcessPyrsRec(/) 

1 addSimplex^i,...,^) 1 for H e Jf < (C i _i,x,-) 

2 for i 4- d + 1 to n 2 do key <- {x ( } U ({jq, . . . ,jc;_i} HH) 

3 do 3 TotalPyrDec (key) 

4 ProcessPyrsRec(z') 



Table 2. Total pyramid decomposition 



first realizes the building of H(x\, . . . ,x n ) (represented by its full dimensional members) 
and the second takes care of the recursion that defines IV a (xi, . . . ,x n ). 

Pyramid decomposition has the virtue of requiring very little memory since one needs 
not store the triangulation in order to produce all the simplicial cones in it. However, there 
is a severe drawback: as above, one must compute the support hyperplanes in Jt?(P) for 
all pyramids encountered. In a "pure" approach, one computes the support hyperplanes 
of the simplicial cones at the bottom of the pyramid decomposition; this is essentially the 
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inversion of the matrix of its generators (see equation (14.11) ). Then one passes them back 
from a pyramid to its "mother", discarding those that fail to have all generators of the 
"mother" in its positive halfspace or have been found previously. These two conditions 
are easily tested. Suppose P is the pyramid to which TotalPyrDec is applied in Pro- 
CESSPyrsRec and G e Jt?(P). Then G E Jt?{x\ ...,xi)\ Jf(x\, . . . ,x ? _i) if and only if 
the following two conditions are satisfied: 

(i) xjeG+ for j=l,...,i-l; 

(ii) Xj G G > for all j = 1 1 such that xj P. 

One should note that pyramids effectively reduce the dimension: the complexity of 
cone(F,;c„) is completely determined by the facet F, which has dimension d— 1. 

While being extremely memory efficient, total pyramid decomposition is usually much 
slower than building the lexicographic triangulation directly. For one of our standard 
test examples (4x4x3 contingency tables, dimension 30 with 48 extreme rays; see O) 
the lexicographic triangulation with respect to the order of generators in the input file has 
2, 654, 272 full dimensional cones. In serial computation on an Intel i7 2600 PC, LexTri- 
ANGULATION computes it in approximately 2 minutes, whereas TotalPyrDec needs 
about 1 1 minutes. The current implementation, described in the next section, reduces the 
serial computation time to 13 seconds. 

Remark 5. Pyramid decomposition is not only extremely useful for the computation of 
triangulations, but also helps in finding support hyperplanes. For them the critical com- 
plexity parameter is | J^ 9< (Q_i,X()| • | > (Q_i,Xj)|, and as in its use for triangulation, 
pyramid decomposition lets us replace a very large product of the sizes of two "global" 
lists by a sum of small "local" products-the price to be paid is the computational waste 
invested for the support hyperplanes of the pyramids that are forgotten later on. Never- 
theless pyramid decomposition leads to a substantial reduction in computing time also 
for support hyperplanes, and Normaliz uses this effect. We illustrate this by computa- 
tion times for the 5x5x3 contingency tables (dimension 55 with 75 extreme rays; see 
0). The cone has 306,955 support hyperplanes. On a Sun xFire 4450 we measured a 
serial computation time of 16,822 seconds if only FindNewHyp is used. The current 
implementation reduces this to 4, 334 seconds. 

3. The current implementation 

Since version 2.7 (and partly since 2.5) Normaliz has combined lexicographic trian- 
gulation with pyramid decomposition. The support hyperplanes and the triangulation are 
extended from one generator to the next only until certain bounds are exceeded. From 
that point on, the algorithm BuildCone described in Table [3] switches to pyramid de- 
composition, and the same mixed strategy is then applied to the pyramids. 

We now use two types of passage to pyramids, a recursive one via ProcessPyrsRec 
and a nonrecursive one via ProcessPyrs. The main reason for this split approach is that 
on the one hand recursion limits the effect of parallelization (as it does in Normaliz 2.7), 
and, on the other hand, the recursive approach nevertheless saves time in the computation 
of support hyperplanes for the top cone. 
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BuildCone (xi, . . . ,x n ; recursion) ProcessPyrsRec(/) 

1 ADDS IMPLEX (xi,...,X d ) 1 for H G M' < (Q- i,X[) 

2 for i 4- d+ 1 to n 2 do key 4- {xj} U ({x\, . . . DH) 

3 do 3 BuildCone (key, true) 

4 if MakePyramidsSupp & recursion 

5 then ProcessPyrsRec(/) PROCESsPYRs(z',/eve/) 

6 else if MakePyramidsTri 1 for H 6 < (Q_i,x,-) 

7 then ProcessPyrs (/, /eve/) 2 do fcey <- {x,} U ({xi, . . . ,x ; -i} n#) 

8 else ExtendTri(z') 3 StorePyr (fcey, /eve/ + 1) 

9 FindNewHyp(z') 

10 if TopCone 

1 1 then EvaluatePyrs (0) 

Table 3. Combining lexicographic triangulation and pyramid decomposition 



The boolean recursion indicates whether the recursive passage to pyramids is allowed. 
For the top cone BuildCone is called with recursion = true. The boolean MakePyra- 
midsSupp combines two conditions: 

(1) while set to false initially, it remains true once the branch ProcessPyrsRec has 
been taken the first time; 

(2) it is set true if the complexity parameter | J^ < (Q_i,x ! -)| • | J^ > (Q_i,x ! -)| exceeds 
a threshold. 

In the nonrecursive passage to pyramids we cut the umbilical cord between a pyramid 
and its mother and just store the pyramid for later evaluation. The nonrecursive call is 
controlled by the boolean MakePyramidsTri that combines three conditions: 

(1) while set to false initially, it remains true once the branch ProcessPyrs has 
been taken the first time; 

(2) it is set true if the complexity parameter 1 34? < (Q_i,x/) | • |E| exceeds a threshold; 

(3) it is set true if the memory protection threshold is exceeded; 

The last point needs to be explained. BuildCone is not only called for the processing 
of the top cone C, but also for the parallelized processing of the stored pyramids. Since 
each of the "parallel" pyramids produces simplicial cones, the buffer in which the sim- 
plicial cones are collected for evaluation, may be severely overrun without condition (3), 
especially if | < (Cj_i,jc,-) | is small, and therefore condition (2) is reached only for large 
|A(xi,. . . ,Xj_i)|. The variable level indicates the generation of the pyramid; for the top 
cone it has value —1, and increases by 1 with each new generation. 

At the end of BuildCone for the top cone C we start the evaluation of the stored 
pyramids as described in Table 0] 

Remark 6. (a) For efficiency Normaliz completely avoids nested parallelization. There- 
fore the parallelization in FindNewHyp and ExtendTri is switched off when the par- 
allelization in EvaluatePyrs is active. On the other hand, these are active when the top 
cone or recursively built pyramids are being processed. 



THE POWER OF PYRAMID DECOMPOSITION IN NORMALIZ 



9 



EvaluatePyrs (level) 

1 if PyramidList[level] = 

2 then return 

3 parallel for P e PyramidList[level] 

4 do 

5 BUILDCONE(P,/a/se) 

6 EvaluatePyrs (level + 1 ) 

Table 4. Evaluation of pyramids 

(b) Despite of considerable efforts we have not found a completely satisfactory solution 
in which pyramids could always be processed recursively and simultaneously in parallel. 
Because of (a) we can only parallelize the pyramids directly produced from the top cone in 
the recursive approach, and then parallelization may be limited by an insufficient number 
of hyperplanes in M >< [Ci-\ y Xi) or, more often, by enormous differences in the sizes of the 
pyramids, so that one of them may be running solo for a long time — recognizing the size 
in advance has turned out difficult. Parallelization in FindNewHyp and ExtendTri is 
then the better solution. 

Moreover, a large pyramid together with its children may produce a huge number of 
simplicial cones and overrun the evaluation buffer. Serial loops can be interrupted at any 
time, and therefore the memory problem cannot arise. 

(c) As soon as BuildCone switches to pyramids, the triangulation A(xi, . . . is 
no longer needed for further extension. Therefore it is shipped to the evaluation buffer. 
The buffer is emptied whenever it has exceeded its preset size and program flow allows 
its parallelized evaluation. (Because of (a) this is not always possible.) 

(d) The strategy for the evaluation of pyramids is similar. If the buffer for level + 1 is 
exceeded, evaluation on that level will be started as soon as possible. Usually this results 
in a tree of evaluations over several levels. 

We add a few minor details of the implementation. 

Remark 7. (a) For nonrecursive pyramids the support hyperplanes arising from the last 
generator need not be computed since they are irrelevant for triangulation and pyramid 
decomposition. 

(b) Simplicial facets of Q_i produce exactly one simplicial pyramid in C;. They are 
treated directly by AddSimplex. 

(c) If the extreme rays of C have been singled out from the given generators x\,...,x n 
before BuildCone is called, then only the extreme rays are used in the pyramid decom- 
position and the lexicographic triangulation. 

(d) If a grading is defined explicitly (see Section S]), then Normaliz orders the generators 
by degree and those of the same degree by input order before building the cone C. This is 
an attempt to cover as much ground as possible by using generators of small degree. On 
the whole, we have reached good results with this choice. 

Remark 8. (Partial triangulation) The idea of pyramid decomposition was born when 
the authors observed that the computation of Hilbert bases usually does not need a full 
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triangulation of C. If a simplicial cone o cannot contribute new candidates for the Hilbert 
basis of C, it need not be evaluated, and if a pyramid consists only of such simplicial 
cones, it need not be triangulated at all. This is the case if ht#(;cj) = 1 in ProcessPyrs. 

The resulting strategy has sometimes striking results and was already described in [|5]. 
It is mentioned here only for completeness. If a full triangulation is not required, then 
ProcessPyrs discards all pyramids of height 1 from further processing. (However, their 
support hyperplanes must be computed if processed recursively.) If followed strictly, the 
recursion will not stop before the simplicial cones at the bottom of the pyramid decompo- 
sition. As for full triangulations, this is usually not optimal. Normaliz therefore switches 
to ExtendTri for pyramids of height > 2 from a certain level on. 



The fast computation of triangulations via pyramid decomposition must be accompa- 
nied by an efficient evaluation of the simplicial cones in the triangulation, which is almost 
always the more time consuming step. 

Let a be a simplicial cone generated by the linearly independent vectors vi,...,v^. 
The evaluation is based on the generator matrix G a whose rows are vi, . . . ,Vd- Before 
we outline the evaluation procedure, let us substantiate the remark made in Section |2]that 
finding the support hyperplanes amounts to the inversion of G a . Let H be the support 

hyperplane of o opposite to vu given by the linear form A; = aue\ -\ Vo-di^d with 

coprime integer coefficients aj. Then 



Thus the matrix (ay) is up to scaling of its columns. Usually the inverse is computed 
only for the first simplicial cone in every pyramid since its support hyperplanes are really 
needed. But matrix inversion is rather expensive, and Normaliz goes to great pains to 
avoid it. 

Normaliz computes sets of vectors, primarily Hilbert bases, but also measures, for ex- 
ample the volumes of rational poly topes. A poly tope P arises from a cone C by cutting C 
with a hyperplane, and for Normaliz such hyperplanes are defined by gradings: a grading 
is a linear form deg : 7h d — > Z (extended naturally to W 1 ) with the following properties: (i) 
deg(x) > for all x e C, x ^ 0, and (ii) deg(Z rf ) = Z. The first condition guarantees that 
the intersection P = CC\A\ for the affine hyperplane 



is compact, and therefore a rational polytope. The second condition is harmless for in- 
tegral linear forms since it can be achieved by extracting the greatest common divisor of 
the coefficients of deg with respect to the dual basis. 

The grading deg can be specified explicitly by the user or chosen implicitly by Nor- 
maliz. The implicit choice makes only sense if there is a natural grading, namely one 
under which the extreme integral generators of C all have the same degree. (If it exists, it 
is of course uniquely determined.) 



4. Evaluation of simplicial cones 



(4.1) 




Ai = {xeR d :deg(x) = \} 
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At present, Normaliz evaluates the simplicial cones a in the triangulation of C for the 
computation of the following data: 

(HB) the Hilbert basis of C, 

(LP) the lattice points in the rational poly tope P = CC\A\, 

(Vol) the normalized volume vol(P) of the rational poly tope P (also called the multi- 
plicity of C), 

(Ehr) the Hilbert or Ehrhart function H(C,k) = \kPC\Z d \,ke Z+. 

Task (Vol) is the easiest, and Normaliz computes vol(i') by summing the volumes 
vol(anAi) where o runs over the simplicial cones in the triangulation. With the no- 
tation introduced above, one has 

vol CTDAi = , ' , , r - 

deg(vi)---deg(v d ) 

For the justification of this formula note that the simplex o DA\ is spanned by the vectors 
v,-/ deg(vf), i = 1, . . . ,d, and that the vertex of the J-simplex 8 = conv(0, oC\A\) has 
(lattice) height 1 over the opposite facet o P\A\ of 8 so that vol(<7 fl A\) = vol(5). 

The remaining tasks depend on the set E of lattice points in the semi-open parallelotope 

par(vi,...,v rf ) = {^ivi H \-q d v d : <q t < 1}. 

For the efficiency of the evaluation it is important to generate E = Z d fl par(vi , . . . , v^) as 
fast as possible. The basic observation is that E is a set of representatives of the group 
Z d /U a where the subgroup U a is spanned by vi, . . . , v^. Thus one finds E in two steps: 

(Rep) find a representative of every residue class, and 

(Mod) reduce its coefficients with respect to the Q-basis vi, . . . , v d modulo 1. 

The first idea for (Rep) that comes to mind (and used in the first version of Normaliz) 
is to decompose Z d /U a into a direct sum of cyclic subgroups Zw,-, i = 1,...,J where 
«1, . . . , Ud is a Z-basis of Z d and _ denotes the residue class modulo U a . The elementary 
divisor theorem guarantees the existence of such a decomposition, and finding it amounts 
to a diagonalization of G a over Z. But diagonalization is even more expensive than matrix 
inversion, and therefore it is very helpful that a filtration ofZ d /U a with cyclic quotients 
is sufficient. Such a filtration can be based on trigonalization: 

Proposition 9. With the notation introduced, let e\, . . . , e d denote the unit vectors in Z rf 
and let X e GL(J, Z) such that XG a is an upper triangular matrix D with diagonal ele- 
ments a i , . . . , ad > 1. Then the vectors 

(4.2) b\e\A Yh d e d , < b t < a u i = 1, • • • , d, 

represent the residue classes in U* /U a . 

Proof. Note that the rows of XG G are a Z-basis of U a - Since \Z d /U(j\ = \ detG a \ = 
a\ • ■ -ad it is enough to show that the elements listed represent pairwise different residue 
classes. Let p be the largest index such that a p > 1. Note that a p is the order of the cyclic 
group Zip, and that we obtain a Z-basis of U' a = U a + Ze p if we replace the p-th row of 

XG a by e p . If two vectors b\e\-\ h b p e p and b\e\ H h b' p e p in our list represent the 

same residue class modulo U a , then even more so modulo U' G . It follows that bi = b\ for 
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i = 1, . . . ,p — 1, and taking the difference of the two vectors, we conclude that b p = b' p as 
well. □ 

The first linear algebra step that comes up is therefore the trigonalization 

(4.3) XG a =D. 

Let be the transpose of G a . For (Mod) it is essentially enough to reduce those ei 
modulo 1 that appear with a coefficient > in (14.21) . and thus we must solve the simulta- 
neous linear systems 

(4.4) G% Xi = e u ai>l, 

where we consider xi and e,- as column vectors. In a crude approach one would simply 
invert the matrix G£ (or G CT ), but in general the number of i such that a, > 1 is small 
compared to d (especially if d is large), and it is much better to solve a linear system with 
the specific multiple right hand side given by (14.41 ) . The linear algebra is of course done 
over Z, using a\-- -ad as a common denominator. Then Normaliz tries to produce the 
residue classes and to reduce them modulo 1 (or, over Z, modulo a\ ■ ■ -ad) as efficiently 
as possible. 

For task (LP) one extracts the vectors of degree 1 from E, and the degree 1 vectors 
collected from all o from the set of lattice points in P = CnA\. For (HB) one first 
reduces the elements of E U {vi , . . . , v^} to a Hilbert basis of o, collects these and then 
applies "global" reduction in C. This procedure has been described in [6|, and nothing 
essential has been added meanwhile. 

The most difficult and mathematically most interesting task is (Ehr). For its solution 
one uses the well-known fact that the Hilbert or Ehrhart series, the generating function 

oo 

H c (t) = Y,H(C,k)t k , 

is a rational function of t. For o one has 
h + h l t + --- + h s t s 

H a[t)= n_ tgl y.j l _ t g d y g; = degv ; -, hj = \{x G E : degx = j}\. 

This follows immediately from the disjoint decomposition 

(4.5) Z d no=\Jx+M a 

xeE 

where M a is the (free) monoid generated by vi, . . . , v^. 

However, one cannot compute Hc(t) by simply adding these functions since points in 
the intersections of the simplicial cones a would be counted several times. Fortunately, 
the intricate inclusion-exclusion problem can be avoided since there exist disjoint decom- 
positions of C by semi-open simplicial cones o\S where S is a union of facets (and not 
just arbitrary faces!) of o. The series H a \ s (t) is as easy to compute as H a (t) itself. Let 
x E E, x = Y,1i v i- Then we define e(x) as the sum of all v, such that (i) qt = and (ii) the 
facet opposite to v ; belongs to S. Then 

(4 - 6) = ( r?«.)... ( i-,«o - 
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This follows from the fact that (x + M a ) \S = e(x) +x + M a , and so we just sum over the 
disjoint decomposition of Z n (<7 \ S) induced by (|4.5I) . (Also see 16, Lemma 1 1].) 

The existence of a disjoint decomposition of C into sets of type a \ S was shown by 
Stanley lf2Tl using the existence of a line shelling of C proved by Bruggesser and Mani. 
Instead of finding a shelling order for the lexicographic tri angulation (which is in princi- 
ple possible), Normaliz 2.0-2.5 used a line shelling for the decomposition, as discussed 
in |[6l. This approach works well for cones of moderate size, but has a major drawback: 
finding the sets S requires searching over the shelling order, and in particular the whole tri- 
angulation must be stored. Koppe and Verdoolaege [18 ] proved a much simpler principle 
for the disjoint decomposition (already implemented in Normaliz 2.7). As a consequence, 
each simplicial cone in the triangulation can be treated in complete independence from the 
others, and can therefore be forgotten once it has been evaluated (unless the user insists 
on seeing the triangulation): 

Lemma 10. Let Oq be a vector in the interior of C such that Oq is not contained in 
a support hyperplane of any simplicial o in a triangulation of C. For o choose as 
the union of the support hyperplanes J^f" c (o,Oc)- Then the semi-open simplicial cones 
O \ S a form a disjoint decomposition ofC. 

See lfT8l for a proof. It is of course not possible to choose an order vector Oq that 
avoids all hyperplanes in advance, but this is not a real problem. Normaliz chooses Oq 
in the interior of the first simplicial cone, and works with a lexicographic infinitesimal 
perturbation 0' c . (This trick is known as "simulation of simplicity" in computational 
geometry; see IfTTTn . If O c G H< (or O c G H>), then 0' c G H< (or 0' c G H>). In the 
critical case Oc G H, we take the linear form A representing H and look up its coordinates 
in the dual basis e*, . . . , e* d . If the first nonzero coordinate is negative, then 0' c G H < , and 
elsetf^G/P. 

At first it seems that one must compute the support hyperplanes of a in order to apply 
Lemma [TOl However, it is much better to solve the system 

(4.7) G%I a = O c . 

The solution I a is called the indicator of a. One has O c G H< (or O c G //>) if If < 
(or If > 0) for the generator v ; - opposite to H (A vanishes on H). Let us call o generic if 
all entries of I a are nonzero. 

If If = — this happens rarely, and extremely rarely for more than one index i — then 
we are forced to compute the linear form representing the support hyperplane opposite of 
v,. In view of (14.11) this amounts to solving the systems 

(4.8) G a x = e u If = 0, 

simultaneously for the lexicographic decision. 

If a is unimodular, in other words, if | detG a | = 1, then the only system to be solved is 
(14.71) . provided that o is generic. Normaliz tries to take advantage of this fact by guessing 
whether o is unimodular, testing two necessary conditions: 

(PU1) Every o (except the first) is inserted into the triangulation with a certain generator 
Xj. Let H be the facet of o opposite to x,-. If ht# (x;) > 1, then o is nonunimodular. 
(The number ht# (x ( -) has been computed in the course of the triangulation.) 
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(PU2) If gcd(deg v\ , . . . , deg v^) > 1, then a is not unimodular. 

If a passes both tests, we call it potentially unimodular. (Data on the efficiency of this 
test will be given in Remark fT27 e). 

After these preparations we can describe the order in which Normaliz treats the trigo- 
nalization ([4731) and the linear systems (1474T ). (14771) and (14/81) : 

(LI) If a is potentially unimodular, then (14.71 ) is solved first. It can now be decided 
whether a is indeed unimodular. 

(L2) If o is not unimodular, then the trigonalization (14.31 ) is carried out next. In the 
potentially unimodular, but nongeneric case, the trigonalization is part of the so- 
lution of (14.81) (with multiple right hand side). 

(L3) In the nonunimodular case, we now solve the system (14.41) (with multiple right 



(L4) If o is not potentially unimodular and not generic, it remains to solve the system 
(14.81) (with multiple right hand side). 

As the reader may check, it is never necessary to perform all 4 steps. In the unimodular 
case, (LI) must be done, and additionally (L2) if o is nongeneric. If o is not even poten- 
tially unimodular, (L2) and (L3) must be done, and additionally (L4) if it is nongeneric. In 
the potentially unimodular, but nonunimodular case, (LI), (L2) and (L3) must be carried 
out. 

Remark 11. (a) If one stores the transformation matrix X of (14.31) and its inverse (for ex- 
ample as a sequence of row exchanges and elementary transformations), then one can 
solve the remaining systems without further trigonalization. However, in general the 
bookkeeping needs more time than it saves as tests have shown. 

(b) The simplicial cones stored in the evaluation buffer are processed in parallel, and 
parallelization works extremely well for them. 

(c) Simplicial cones of height 1 need not be evaluated for (HB) and (LP); see Remark[8] 

We conclude this section with a brief discussion of the computation and the represen- 
tation of the Hilbert series by Normaliz. The reader can find the necessary background in 
JH Chapter 6]. 

Adding the Hilbert series (14.61) is very simple if they all have the same denominator, 
for example in the case in which the generators of C (or at least the extreme integral 
generators) have degree 1. For efficiency, Normaliz first forms "denominator classes" in 
which the Hilbert series with the same denominator are accumulated. At the end, the class 
sums are added over a common denominator that is extended whenever necessary. This 
yields a "raw" form of the Hilbert series of type 



hand side). 



(4.9) 



Hc(t) 



R(t) 



R(t)eZ[t], 



(l-Fi 



)...(i-^r 



whose denominator in general has > d factors. 

In order to find a presentation with d factors, Normaliz proceeds as follows. First it 
reduces the fraction to lowest terms by factoring the denominator of (14.91 ) into a product 



THE POWER OF PYRAMID DECOMPOSITION IN NORMALIZ 



15 



of cyclotomic polynomials: 

(4.10) H c (t) = Z{t) , Z(t)eZ[t], &j\Z(t), 

which is of course the most economical way for representing Hc(t) (as a single fraction). 
The orders and the multiplicities of the cyclotomic polynomials can easily be bounded 
since all denominators in (14.61) divide (1 — t ) d where £ is the least common multiple of 
the degrees degx, . So we can find a representation 

Fit) 

(4-U) H c(t) = - / r^4; r, F(t)eZ[t], 

LW (1 — ) - - - (1 -t e d) w u 

in which is the least common multiple of the orders of the cyclotomic polynomials that 
appear in (14.101) . ej-i is the least common multiple of the orders that have multiplicity 
> 2 etc. Normaliz produces the presentation (14.1 II ) whenever the degree of the numerator 
remains of reasonable size. 

It is well-known that the Hilbert function itself is a quasipolynomial: 

(4.12) H(C,k)=q (k)+q l (k)k+---+q d _ l (k)k d - 1 , k>0, 

where the coefficients qj(k) E Q are periodic functions of k whose common period is the 
least common multiple of the orders of the cyclotomic polynomials in the denominator of 
(14.101) . Normaliz computes the quasipolynomial, with the proviso that its period is not too 
large. It is not hard to see that the periods of the individual coefficients are related to the 
representation (14.1 11 ) in the following way: is the common period of the coefficients 
qa-i,.--,qd-k' The leading coefficient q d -\ is actually constant (hence e\ = 1), and 
related to the multiplicity by the equation 

volfP) 

(4.13) q d _! 



(d-l] 



Since q d -\ and vol(P) are computed completely independently from each other, equation 
(14.131) can be regarded as a test of correctness for both numbers. 

The choice (14.1 II) for Hc(t) is motivated by the desire to find a standardized represen- 
tation whose denominator conveys useful information. The reader should note that this 
form is not always the expected one. For example, for C = with deg(ei) = 2 and 
deg(e2) = 3, the three representations (I4.91 l- (|4.1 II ) are 



1 1 l-t+t 



2 



(i-' 2 )0-' 3 ) d 2 C 2 C 3 (i-0(i-^ 6 )' 

Actually, it is unclear what the most natural standardized representation of the Hilbert 
series as a fraction of two polynomials should look like, unless the denominator is (1 — t) d . 
Perhaps the most satisfactory representation should use a denominator (1 — t Px ) ■ ■ ■ (1 — 
t Pd ) in which the exponents p\ are the degrees of a homogeneous system of parameters 
(for the monoid algebra K\L d DC] over an infinite field K). At present Normaliz cannot 
find such a representation (except the one with the trivial denominator (1 —t e ) d )), but 
future versions may contain it. 
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5. Computational results 

A driving force for the recent improvements in Normaliz that we have described in the 
previous sections was the desire to compute the volumes and Ehrhart series of certain 
polytopes related to combinatorial voting theory. 

5.1. Voting schemes and volumes of rational polytopes. We briefly sketch the connec- 
tion between rational polytopes and combinatorial voting theory, referring the reader to 
US, Il20ll or Il23l for a more extensive treatment. We consider the three voting schemes 
discussed in lEOll . 

Consider an election in which each of the k voters fixes a linear preference order of n 
candidates. In other words, voter i chooses a linear order j\ y t ■ ■ ■ ■ H j n of the candidates 
1 , . . . , n. Set N = n\. Counting the preference orders gives an AMuple (vi , . . . , v^) in which 

v p is the number of voters that have chosen the preference order p. Then v\-\ \-vn = k, 

and (vi, . . . , vn) can be considered as a lattice point in the positive orthant of W^, more 
precisely, as a lattice point in the simplex 

% {n) = DA k = DAi) = 

where is the hyperplane defined by x\ H \-xn = k, and = is the unit 

simplex of dimension N — I naturally embedded in A^-space. All further discussion is 
based on the Impartial Anonymous Culture assumption that all lattice points in the simplex 
^4 have equal probability of being the outcome of the election. 

We fix a specific outcome v = (vi, . . . ,vjv). Let us say that candidate j beats candidate 
/ with respect to v if 

(5.1) \{i:jy i j':i=l,...,k}\>\{i:fy i j:i=l,...,k}\. 

As the Marquis de Condorcet observed, the relation "beats" is nontransitive in general, 
and one must ask for the probability of Condorcet's paradoxon, namely an outcome with- 
out a Condorcet winner where candidate j is a Condorcet winner if j beats all other 

candidates /. Let C^' (j) denote the probability that candidate j is the Condorcet winner, 

(n) 

and C k the probability that there is a Condorcet winner. By symmetry and by mutual 

exclusion cj^ = nCi (1). Usually the number k of voters is very large, and therefore one 
is mainly interested in the limit 



CW = lim C| n) = n lim d n) (l) = nd n )(l). 

Let us fix candidate 1 . It is not hard to see that the n — 1 inequalities (15.11) for j = 1 
and / = 2, .. . ,n constitute homogeneous linear inequalities in the variables vi, . . . ,vjv. 

Together with the inequalities Vj > they define a semi-open subpolytope ?t"' of ^jf 1 ' '. 
Then 

, |« (B) nZ^| vol^, W - M 
(5.2) C«(l) = lim i-^ 1 = 7-7 = vol^ W 



k- 



\% (n, r\Z N \ vo\% [ 



where denotes closure and = . For the validity of (15.21 ) note that we work with 
the lattice normalized volume in which the unit simplex has volume 1 . 
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In the case of two candidates Concordet's paradox cannot occur (if one excludes draws), 
and for 3 candidates the relevant volume is not hard to compute. The situation changes 
significantly for 4 candidates since c ta^ has dimension 23 and 234 vertices. As a subpoly- 
tope of ^( 4 ), ^( 4 ) is cut out by the inequalities A ; (v) > 0, i = 1, 2, 3 whose coefficients are 
in the first 3 rows displayed in Table[5] For the assignment of indices the preference orders 
are listed lexicographically, starting with 1 y 2 y 3 y 4 and ending with 4 y 3 y 2 y 1 . 
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Table 5. Inequalities for ^ and 



The lexicographic triangulation used by Normaliz has (only) 1,473,107 simplicial 
cones. Normaliz computes 

1717 



vo 



8192 

in a few seconds. It follows that C< 4 ) = 1717/2048 « 0.8384. According to [H2, this 
value was first determined by Gehrlein [[141 . 

The simplest way out of the dilemma that there may not exist a Condorcet winner 
is plurality voting: candidate j is the plurality winner if j has more first places in the 
preference orders of the voters than any of the other n — 1 candidates. The Condorcet 
efficiency of plurality voting (and similarly of other voting schemes) is the conditional 
probability that the Condorcet winner is elected, provided it exists, as k — > °°. Therefore 
one must compute the probability of the event that candidate j is the Condorcet winner, 
but candidate / 7^ j is the plurality winner. By symmetry, one can assume 7 = 1 and 
ji _ 2 rp^ £ sem j_ p en polytope whose lattice points represent this unexpected 

An) 

outcome is cut out from ^ ' by n — 1 further inequalities saying that 2 has more first 
places than the other n — 1 candidates. Thus one obtains 

C(») - n (n- 1) vol &( n ) 

as the Condorcet efficiency of plurality voting where ^( n ) = ^ >< f\ 

The extra 3 inequalities A,-(v) > 0, i = 4, 5, 6, given in the last 3 lines of Table [6] increase 
the complexity of the polytope ^( 4 ) enormously. It has 3928 vertices, and the triangu- 
lation increases to 347,225,775,338 simplicial cones. Nevertheless, Normaliz computes 
the volume and the Ehrhart series in acceptable time. We have obtained 

3694037 1 852901 6355068 149 1 



vo l £*(4) 



2054269543278 1 8240000000000 
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so that the Condorcet efficiency of plurality voting turns out to be 

d 4 )-12vol^( 4 ) 10658098255011916449318509 nriAnr 
CW ~ 14352135440302080000000000 ~ ' 

in perfect accordance with EOl . 

Additional support (and tests for Normaliz) can be added by a chain of volume compu- 
tations as follows. We start from the unit simplex Pq = and set 

P+ = {xefy {4) : Ai(jc),...,Af(jc) >0} and Pf = P+ 1 n{x: A;(x) <0}, i=l,...,6. 

Then P^ = and we must have 



voli^Ij = vol/f + vol/T, i= 1,...,6. 

The volumes computed by Normaliz satisfy all these equations, as they should. We start 
from i = 3 since P^ = 



volP+ 



1717 
8192' 



, , 418988423262545 , 2982999236660911 

vol Pf = , vol P7 



16231265527136256' 4 16231265527136256' 

, , 1622886339180775733501803 , 365671997787943700091947 
volP = volP = 

5 77035107872931840000000000' 5 77035107872931840000000000' 

, , 3694037185290163550681491 , 1900979157575715215969951 
volP = volP - = 

6 205426954327818240000000000' 6 616280862983454720000000000' 

The largest triangulation of 463,613, 250,401 simplicial cones was produced by PjT. 

The last problem discussed in [|20l is plurality voting versus plurality cutoff. It works 
as follows. In the first round of the election the two top candidates in plurality voting 
are selected, and in the second round the preference orders are restricted to these two 
candidates. In order to model this situation by inequalities one must fix an outcome of the 
first round that involves all n candidates, say 1, . . . ,n in this order. This condition gives 
rise to n — 1 inequalities. Then the n-th inequality expresses that 2 is the winner of the 
second round, despite the fact that 1 was the winner of the first round. The volume of 
the corresponding polytope gives the probability of this event. By mutual exclusion and 
symmetry, we must multiply the volume by n! in order to obtain the probability for the 
event that the winner of the first plurality round looses after cutoff. 

As a subpolytope of the polytope is defined by the inequalities in Tabled 



1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 
1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 
1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 

-1-1-1-1-1-1 1 1 1 1 1 1-1 1-1 1-1 1-1 1-1 1-1 1 



Table 6. Inequalities for 
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It has 1 872 vertices and the triangulation computed by Normaliz has 257, 744, 341 , 008 
simplicial cones. The volume is 

, „ U ) 2988379676768359 

vol £2 y 1 



292162779488452608 

The total probability of the failure of the winner of the first round is therefore 24 • vol i?( 4 ) 
~ 0.2455, in accordance with the results of [|2~0l for this model. Further support has 
been given by De Loera, Dutra, Koppe, Moreinis, Pinto and Wu in lfl2l . where LattE 
Integrale IITSI was used for the volume computation. (Additionally we have applied the 
same verification as for &^). 

5.2. Ehrhart series and quasipolynomials. Normaliz has not only computed the vol- 
umes, but also the Ehrhart series and quasipolynomials for the closures of the semi-open 

polytopes <t?( 4 \ ^ and The Hilbert series of is the rational function with 

numerator 

1 + 1 + 5t 1 + 133? 2 + 363? 3 + 458 1? 4 

+ 8655? 5 + 69821? 6 + 100915? 7 + 596834? 8 + 697232f 9 

+ 3255226? 10 + 3176870? 11 + 12235441f 12 + 10182887? 13 + 33268048? 14 

+ 23917200? 15 + 67509138? 16 + 42243510? 17 + 104272000f 18 + 56990048; 19 

+ 123966919? 20 + 59177761? 21 + 1 13925878; 22 + 47336170? 23 + 80758791? 24 

+ 28993857? 25 + 43770180? 26 + 13415068? 27 + 17837843? 28 + 4580485? 29 

+ 5320122? 30 + 111 1974? 31 + 1 1 13216? 32 + 180850/ 133 + 152 891? 34 

+ 17845? 35 + 12346f 36 + 890f 37 + 481/ 38 + 15? 39 + 6t 40 . 

and denominator 

(i-0(i-* 2 ) 14 (i-* 4 ) 9 

Numerator and denominator are coprime. 

— (4) 

The quasipolynomial of c € has period 4. We give the numerators of its first two and 
its last two coefficients; the denominator is d = 6939597901822221635907747840000: 

q (0) = 6939597901822221635907747840000/J= 1, 

qi (0) = 20899225 148336747959025664000000/ d, 

422(0) = 422(2) = 15982652919/J, 

423(0) = 423(1) = 423(2) = 423(3) = 56262656/J = vol^( 4 )/23!, 

q (l) = 2034750310223351797008092160000/J, 

qi (l) = 7092764342142539187142971648000/J, 
422(1) = ^2(3) = 15528493056/ d, 

q (2) = 6933081849299152199775682560000/J, 

qi (2) = 2089245531 1735756236854919168000/J, 

?o(3) = 2034750310223351797008092160000/J, 
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4i (3) = 7092764342142539187142971648000/d. 

The reader may have noticed that the coefficients of q{\) and q(3) coincide, as far as 
listed. In fact, these two polynomials are equal. 

For the other two polytopes we only list the denominators for the representation (14.1 II ) 
(with non-coprime numerators in both cases): 

W {4) : (i -0(1 -t 2 ) 2 (\ -t 4 ) 5 (i -t n )\\ -t 24 )(i -t m ) l \ 

£ {4) : (l-t)(l-t 2 ) 2 (l-t 4 ) 5 (l-t 12 ) 16 . 

(4) — (4) 

So has period 120 and has period 12. On request the authors will provide full 
sets of data. 

For the exact counting of the frequency of Condorcet's paradoxon, and similarly for 
the other events considered, one must compute the Ehrhart series of the semi-open poly- 
topes. At present Normaliz cannot do this directly so that one depends on cumbersome 
inclusion/exclusion for the semi-open polytopes. We plan an extension for semi-open 
polytopes (or cones) the next version. 

One should note that semi-open polytopes present an inherent difficulty since in general 
a disjoint decomposition into sets o \ S as discussed in Section |4] does not exist. (Recall 
that S is a union of facets, not just arbitrary faces.) Suppose that all extreme integral gen- 
erators of the cone C have degree 1 . Then all the Hilbert series H a \ s have denominator 
( 1 — t) d and a numerator polynomial with nonnegative coefficients. If a disjoint decompo- 
sition into sets o \ S exists, the resulting Hilbert series must have a numerator polynomial 
with nonnegative coefficients as well. However, one can easily find semi-open polytopes 
for which this is not the case: remove two opposite edges from the unit square; then the 
Ehrhart series of the remaining semi-open polytope is {It 2 — t 3 ) / (I — t) 3 . 

5.3. Computation times. Table |7] gives an indication of the computation times to be 
expected for the volumes and Ehrhart series of the rational polytopes discussed in the 
previous subsection. The times have been collected at different stages of the development 
on a SUN xFire 4450 with 20 threads (of the maximal number of 24). However, these 
have differed only little in this respect, except that the computation for has become 
about 20% faster now. 



Polytope 


computation 


triangulation size 


real time 


parllelization 




Ehrhart series 


1,473,107 


00:00:30 h 


serial 




Ehrhart series 


347,225,775,338 


292:50:22 h 


1981% 


<g(4) 


Ehrhart series 


257,744,341,008 


175:11:26 h 


1991% 




volume 


383,986,938,515 


126:06:57 h 


1953% 


<2(4) 


triangulation 


271,164,705,162 


18:50:47 h 


1779% 



Table 7. Computation times 
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Remark 12. (a) The size of the lexicographic triangulation depends on the order in which 
the extreme rays are processed. The polytopes in the table above are defined by their 
support hyperplanes, and therefore Normaliz first computes the extreme rays from them. 
The order used in the computations mentioned in Table |7] is not necessarily identical 
with the order produced by the most recent version. In it we have eliminated any of the 
unpredictable effects of parallelization in the function FindNewHyp (see Section |2]). 
See also Remark Hd). 

That the order of the generators has some influence is shown by the two computations 
ofgW in Table H 

(b) The table shows that the times needed for (i) pure triangulation, (ii) volume com- 
putation and (iii) Ehrhart series are in approximate proportion 1:5:10 for this class of 
cones. 

(c) For the calculation of the Hilbert bases of the cones defined by the polytopes ^( 4 ), 
<^( 4 ) and ^ one should use the dual mode of Normaliz. Then the Hilbert basis calcula- 
tions are a matter of seconds. The Hilbert bases have the following numbers of elements: 
242 for tf (4 \ 25192 for 9621 for 

(d) Normaliz needs relatively little memory. All the computations mentioned run stably 
with < 1 GB of RAM. 

(e) From the Ehrhart series calculation of we have obtained the following statis- 
tics on the types of simplicial cones: 61, 845,707, 957 are unimodular, 108,915, 272, 879 
are not unimodular, but satisfy condition (PU1), of which 62, 602, 898, 779 are potentially 
unimodular. This shows that condition (PU2) that was added at a later stage has a satisfac- 
tory effect. (The number of potentially unimodular, but nonunimodular simplicial cones 
is rather high in this class.) The average value of | detG a | is ~ 10. 

The number of nongeneric simplicial cones is 129,661,342. The total number s of 
linear systems that had to be solved for the computation of the Ehrhart series is bounded 
by 516, 245, 872, 838 < s < 516, 375, 534, 180. 

The total number of pyramids was 80,510,681. It depends on the number of parallel 
threads that are allowed. 

5.4. The exploitation of symmetry. The elegant approach of Schiirmann in [20] for the 
computation of the volumes of c io^\ ^( 4 ) and uses the high degree of symmetries 
of these polytopes. If certain variables v ;i , . . . , vt u occur in all of the linear forms given in 
Tables [5] and [6l then any permutation of them acts as a symmetry on the corresponding 

polytope, and the variables v ;i , . . . , Vi u can be replaced by their sum v ;i H h v; 1( in them. 

(The polytopes have further symmetries.) The substitution can be used for a projection 
into a space of much lower dimension, mapping the polytope P under consideration to a 
polytope Q (this requires that the grading affine hyperplane A\ is mapped onto an affine 
hyperplane by the projection). Instead of counting the lattice points in kP one counts the 
lattice points in kQ weighted with their number of preimage lattice points in kP. This 
amounts to the consideration of a generalized Ehrhart function 



k^ £ f(x). 

xekQnZ d 
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The theory of generalized Ehrhart functions has recently been developed in several papers; 
see (191 , [QQ|, • An extension of Normaliz to the computation of generalized Ehrhart 
functions and their generating functions is envisaged. 

In (201 . only the leading term of the polynomial / is used. Integration with respect to 
Lesbesgue measure then yields the volume. 

5.5. Previous challenging computations by Normaliz. We conclude by listing some 
more performance data of previous computations in Table [8] For the 5x5x3 contin- 
gency tables see (51, and the computations of the statistical rank models were done for 
(221 . These computations, for which the triangulations are very large, but not as mon- 
strous as those of <^( 4 ) and £?( 4 \ are now doable in comfortable time. The Hilbert basis 
computations show the efficiency of partial triangulations (see Remark [8]). The compu- 
tations were done on our SUN xFire 4450 with 20 parallel threads at various stages of 
the development. The degree of parallelization varies somewhat because of the problems 
discussed in Remark |6£b). The computation times for triangulation, volume and Hilbert 





dim 


extr 


computation 


triangulation 


real time 


5x5x3 

contingency tables 


43 


75 


Hilbert series 
Hilbert basis 


9,248,527,905 
448,64 


07:07:30 h 
00:20:31 h 


linear rank 
model for S(, 


16 


720 


Hilbert series 
Hilbert basis 


5,745,903,354 
7,783,191 


02:23:04 h 
00:09:56 h 


ascending rank 
model for S5 


27 


120 


Hilbert series 
Hilbert basis 


20,853,141,970 



07:28:11 h 
00:00:01 h 



Table 8. Previous challenging computations 



series of the 5x5x3 contingency tables are in approximate proportion 2:5:6. This is 
not surprising since almost all simplicial cones are unimodular. 
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